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1  Introduction 

The  primary  topic  to  be  discussed  in  this  paper  is  the  direct  numerical  simulation 
of  those  highly  compressible  turbulent  flows  which  can  be  described  by  the 
single-fluid  Navier-Stokes  equations  within  the  constraint  of  periodic  boundary 
conditions.  While  modern  theoretical  analyses  of  compressible  turbulent  flows 
began  with  the  work  of  Moyal  [1],  their  direct  numerical  simulation  is  much 
more  recent.  Following  the  seminal  work  of  Passot  and  Pouquet  [2],  many  other 
authors  have  also  used  Fourier  methods,  for  example,  Erlebacher  tt  ai,  [3], 
Blaisdell  et  al.,  [4],  Sarkar  et  ai,  [5],  Kida  and  Orszag  [6],  and  Zang  ei  al.,  [7]. 

Here,  Fourier  methods  are  extended  in  two  principal  ways.  First,  we  de¬ 
scribe  a  technique  which  has  proven  useful  previously  [8]:  the  governing  equa¬ 
tions  are  cast  in  a  form  where  the  important  physical  variables  are  not  the  fluid 
density  and  temperature  directly,  but  rather  their  natural  logarithms;  this  en¬ 
sures  adherence  to  a  physical  constraint  of  positive-definiteness  which  may  be 
computationally  violated  in  non-logarithmic  formulations,  leading  to  numerical 
instability.  Second,  bulk  viscosity  is  utilized,  both  to  model  polyatomic  g2ises 
more  accurately  and  concomitantly  to  ensure  numerical  stability  in  the  presence 
of  strong  shocks. 

The  efficacy  of  using  logarithmic  variables  and  physical  values  of  bulk  viscos¬ 
ity  will  be  shown  through  several  numerical  examples,  [n  particular,  logarithmic 
variables  will  allow  for  the  simulation  of  supersonic  homogeneous  turbulent  flow, 
while  bulk  viscosity  will  enable  the  resolution  of  shock  structure.  These  numer¬ 
ical  examples  will  be  followed  by  a  conclusion  where  extensions  of  the  current 
work  will  be  discussed. 


2  Basic  Equations 

The  basic  equations  of  compressible  fluid  dynamics  may  be  expressed  as 


dp  _ 

0 

(1) 

dpvi  _ 

-Vp-b  V  •  |pVu-l-  +  IV  •  uj 

(2) 

dp 

jj+V  pu  = 

-(7  -  l)pV  •  u  -t-  (7  -  1)V  •  (kVT) 

+  (7  -  1)  [|  -TijTij  -b  /<b(V  •  uf 

(3) 

Here,  r,j  =  diUj  +  djUi  —  ‘2/36ijV  ■  u  and  I  =  is  the  unit  dyadic. 

The  equation  of  state  will  be  that  of  an  ideal  gas,  p  =  RpT/m,  where 
R/m  =  Cp  —  Cv\  /?  is  the  ideal  gas  constant,  m  is  the  molecular  weight  of 
the  gas  and  Cp  and  c„  are  the  specific  heats  at  constant  pressure  and  volume, 
respectively.  Here,  we  will  consider  polytropic  gases,  i.e.,  gases  such  that  the 
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specific  heats  and  their  ratio  7  =  Cp/c„  are  taken  to  be  constants.  In  this  case, 
the  speed  of  sound  c  satisfies  =  jRT/m  =  jp/p- 

We  can  non-dimensionalize  the  equations  (l)-(3)  in  terms  of  reference  values 
Po,  To,  and  Ug.  Thus,  using  =  yRT„/m  and  po  =  po<^\h  l^he 

following,  dimensionless  quantities  have  a  superscript  attached,  t.g.,  p*): 

P  =  PoP" 

U  =  UqU’ 

T  =  ToT 

P  =  PoP*  -  PocIp^T* /y  (4) 


Now  we  choose  a  length  scale  Lo  so  that  t  =  Igt*  where  to  =  Lo/ug.  Then 
dividing  (1)  by  po/to,  (2)  by  poUg/io,  and  (using  R/m  =  Ct,{y  -  1))  dividing  (3) 
by  Po/to  gives: 


Here,  Mg  =  Ug/Cg  is  the  reference  Mach  number  and  the  dimensionless  transport 
coefficients  are: 


/«*  =  nl{,pgU„Lg) 

PB  =  PB/ipoUglg) 

K*  =  K/{Cgp„UgLg)  (8) 


and  the  t*j  are  the  same  as  previously  defined,  except  now  in  terms  of  the 
corresponding  dimensionless  quantities. 

In  (8),  we  recognize  /<*  =  l/Re,  i.e.,  the  dimensionless  shear  viscosity  is  the 
inver.se  of  the  Reynolds  number  Re.  Although  the  transport  coefficients  which 
appear  in  the  fluid  equations  are  generally  dependent  on  p  or  T  or  both,  here 
we  will  operate  under  the  assumption  that  these  coefficients  are  all  constant. 
Then  the  Prandtl  number  is  Pr  =  Cpfi/n  —  yp*  /k*  and  the  ratio  of  bulk  to 
shear  viscosity  (i  =  pb/p  =  Pb/p"  "'•I'  *1'’°  constant. 

Using  the  quantities  developed  above,  we  can  write  the  dimensionless  fluiil 
equations  in  a  more  compact  form.  First,  we  will  drop  the  from  the  dimen¬ 
sionless  quantities,  remembering  that  we  are  now  dealing  with  non-dimensional 
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equations.  Second,  we  will  choose  the  reference  velocity  as  «„  =  Co  so  that 
Mo  =  1.  Then,  since  7,  //,  fig,  and  k  are  all  constants,  the  equations  become 


^  +  V-^u  =  0  (9) 

dpu  _  I  ^  ^ 

-;^  +  Vpuu  =  — VpT 

01  7 

+  /tV-  Vu+Q  +  /?^IV  u  (10) 

^  +  V  pTu  =  -(j  -  l)pTV  ■u  + kV^T 

+  7(7-1  )/t  ^  Tij  Tij  +  ^(  V  •  u)^  (11) 

(The  dimensionless  pressure  is  p  =  pT /y.)  These  are  the  set  of  bcisic  equations 
with  which  we  will  simulate  the  motion  of  a  compressible  fluid.  In  particular, 
by  judiciously  choosing  appropriate  values  for  /i,  /?,  and  k,  we  will  be  able  to 
simulate  both  turbulence  and  shocks,  as  well  as  their  interaction. 

3  Logarithmic  Variables 

Although  the  variables  p  and  T  may  take  values  only  between  0  and  00,  it 
is  always  possible  in  a  dis.;rete  numerical  simulation  to  inadvertently  assign 
these  variables  non-positive  values.  When  this  occurs,  an  instability  generally 
arises  which  stops  the  simulation.  One  sure  way  to  avoid  this  is  to  express 
the  basic  equations  in  terms  of  logarithmic  variables:  A  =  In  p  and  o-  =  InT. 
This  is  particularly  appropriate  when  we  wish  to  simulate  fluid  flows  in  which 
compressibility  plays  an  important  role,  as  in  the  case  where  shocks  are  present. 
Since  this  is  the  case  at  hand,  we  will  use  these  logarithmic  variables  here. 

Placing  p  =  e^  and  T  =  e”  into  (9)-(Il)  yields  the  basic  non-dimensional 
equations  in  a  logarithmic  formulation: 

-Vu  (12) 

--c'^V(A  +  (7) 

7 

+  pe-^V-  Vu+ Q  +  /^^IV  u  (13) 

-(7-  l)V  u+Ke-'[VV  +  (V(T)2] 

+  7(7  -  + /?(V  u)^  (14) 
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It  will  be  seen  presently  that  these  equations,  though  not  in  conservative  form 
[9],  are  well  suited  to  simulating  highly  compressible  fluid  motion. 


4  Transport  Coefficients 

Along  with  the  value  given  to  the  ratio  of  specific  heats  7,  a  very  important 
feature  of  any  simulation  of  fluid  dynamics  is  the  set  of  values  assigned  to 
the  dimensionless  transport  coefficients  /i,  ftp,  and  k.  Although  all  of  these 
coefficients  vary  with  temperature  and  pressure,  here  they  will  be  assumed  to 
be  constant.  An  investigation  of  the  effects  that  their  dependence  on  density 
and  temperature  may  have  will  be  deferred. 

The  values  assigned  to  these  dimensionless  coefficients  depend  on  both  the 
physical  references  Uo,  Lo,  and  T„,  and  on  numerical  constraints,  such  as  grid 
spacing.  It  is  well  known,  for  example,  that  the  ratio  of  the  largest  dynamically 
important  scales  to  the  smallest  is  proportional  to  [10].  Assuming  that 

the  constant  of  proportionality  is  about  one,  and  since  current  computers  limit 
the  total  number  of  points  on  a  grid  to  about  10®,  we  have  Re  <  500  for  a  three- 
dimensional  grid  and  Rc  <  10'*  for  a  two-dimensional  grid,  approximately.  Note 
that  this  is  the  Reynolds  number  based  on  a  characteristic  physical  length,  and 
not  the  microscale  Reynolds  number,  which  is  based  on  an  average  turbulent 
eddy  size.  Other  numerical  constraints,  such  as  a  limit  on  available  computa¬ 
tional  time,  reduce  grid  sizes  from  their  potential  maximum  and  further  reduce 
the  maximum  Re  which  may  be  considered. 

However,  once  /<  =  1/Re  is  set,  the  other  transport  coefficients  follow  some¬ 
what  directly.  Since  Pr  =  7k//(  a;  1  and  since  1  <  7  <  5/3  for  common  gases, 
we  see  k  %  ft.  The  value  of  ftp,  or  equivalently  (i  =  ftp/ ft,  is  highly  dependent 
on  whether  the  gas  is  polyatomic  or  not;  although  monatomic  gases  have  w  0, 
polyatomic  gases  range  from  /?  =»  1  for  air  to  (i  ss  30  for  molecular  hydrogen  to 
()  «  10^  for  carbon  dioxide  [1 1]. 

5  Numerical  Method 

Here,  we  will  use  a  Fourier  pseudospectral  technique  [12]  to  numerically  solve  the 
logarithmic  variable  equations  (r2)-(  14).  In  this  method,  the  physical  variables 
A,  11,  and  cr  are  expanded  in  terms  of  discrete  Fourier  .series,  c.g., 

A(x)  =  Y. 

lk|<N/2 

The  argument  of  the  variable  will  be  used  to  denote  whether  it  is  in  x-space: 
A(x),  or  in  k-space;  A(k).  In  the  above  equation,  n  is  the  spatial  dimension  and 
N  is  the  number  of  points  on  the  numerical  grid  in  one  dimension  (in  2D  the 
grid  h  N  X  N  and  in  3D  it  \s  N  x  N  x  N).  When  discrete  forward  or  inverse 


4 


Fourier  transforms  are  needed,  they  are  found  using  a  fast  Fourier  transform 
(FFT)  suhroutine. 

Notice  that  the  maximum  value  of  any  component  of  k  is  strictly  less  than 
N l'2.  The  reason  is  that  a  Fourier  coefficient  which  hcus  one  of  the  components 
of  k  equal  to  N /‘I  has  only  a  real  part,  as  far  cis  the  FFT  is  concerned.  Multiply¬ 
ing  this  component  by  rk  yields  the  corresponding  components  of  its  gradient, 
for  example,  which  have  only  imaginary  parts.  An  FFT,  however,  ignores  the 
imaginary  part  of  such  a  component,  so  that  its  contribution  to  any  derivative 
is  always  zero.  It  is  therefore  prudent  to  never  include  such  components  in  any 
FFT-based  method  for  the  numerical  solution  of  differential  equations. 

When  the  Fourier  expansions  of  the  physical  variables  are  placed  into  the 
small  set  of  partial  differential  equations  (FDEs)  (1‘2)-(14),  the  result  is  a  large 
set  of  ordinary  differential  equations  (ODEs),  one  equation  for  every  value  of 
k  which  the  FFT  utilizes.  Although  either  the  set  of  PDEs  in  x-space  or  the 
set  of  ODEs  in  k-space  can  be  time-integrated,  here  it  is  the  ODEs  which  are 
integrated  forward  in  time.  The  reason  is  that  the  numerical  arrays  in  k-space 
contain  fewer  nonzero  elements  than  the  arrays  in  x-space.  Since  |k|  <  N /‘I 
and  since  the  arrays  in  either  space  contain  about  elements,  then  the  ratio 
of  nonzero  array  elements  in  k-space  to  x-space  is  %/A  =  0.79  in  2-D  and 
t/6  =  0.52  in  3D  (i.c.,  the  ratio  of  the  area  of  a  circle  to  the  area  of  the  square 
which  just  encloses  it,  and  the  volume  of  a  sphere  to  the  volume  of  a  cube  which 
Just  encloses  it,  respectively).  Thus,  one  may  always  reduce  the  size  of  a  majority 
of  the  arrays  in  a  k-space-based  computer  code  so  that  they  are  minimal  anti 
yet  contain  all  essential  information.  When  it  is  necessary  to  perform  an  FFT  to 
x-space,  a  minimal  array  is  mapped  into  a  full-sized  FFT  array;  this  is  needed 
only  to  evaluate  products  of  two  x-space  arrays,  so  only  two  full-sized  FFT 
arrays  are  actually  needed.  Since  the  total  number  of  arrays  needed  in  a  typical 
simulation  are  usually  much  greater  than  two,  memory  requirements  in  2-D 
may  be  reduced  by  about  20%  and  in  3D  by  about  50%.  (The  argument  in  this 
paragraph  also  pertains  to  fully  spectral  methods  where  |k|  <  kmar  <  N/2.) 

Although  minimal  arrays  are  critical  in  order  to  maximize  N ,  on  modern 
supercomputers  with  very  large  memories  such  a  large  value  of  N  may  lead  to 
prohibitively  long  run  times.  Since  grid  sizes  are  kept  to  a  reasonable  value 
in  the  work  to  be  presented  in  this  paper,  minimal  arrays  are  not  used  here 
(though  they  have  been  implemented  previously  [13]).  Nonetheless,  this  is  the 
motivation  for  .solving  the  ODEs  in  k-space. 

The  time-integration  method  used  here  was  a  ‘third-order  partially  corrected 
Adams-Bashforth  scheme’  [14].  The  time-step  size  was  variable  and  inversely 
proportional  to  the  largest  absolute  modal  value  (Rmai)  from  the  right-hand- 
side  of  (12)-(14)  at  each  iteration: 


_  f  R-inar  if  Rn 

\  f»  if  R„ 


where  typically  10  ^  <  a  <  10  '.  This  method  of  determining  A/  automatically 


satisfies  the  necessary  stability  conditions  [2]. 

In  the  pseudospectral  method  presented  here,  shock  structure  is  resolved. 
This  is  done  by  considering  cases  where  fi  —  is  large  enough  so  that 

whatever  shocks  occur  are  of  a  naturally  limited  steepness  [15];  for  example, 
the  gas  under  consideration  can  be  Ho  (molecular  hydrogen)  or  some  mixture 
containing  or  the  characteristic  length  Lo  can  be  cussumed  small  for  an 
arbitrary  polyatomic  gas.  In  contrast,  numerical  solutions  of  the  Euler  ecpjations 
require  methods  that  use  shock-capturing  or  shock-fitting  [12,  pp.  255-273]. 

6  Initial  Conditions 

In  order  to  begin  a  simulation,  the  initial  values  of  A,  ii,  and  cr  need  to  be 
specified  (but  not  the  boundary  conditions,  which  are  periodic).  The  initial 
turbulent  velocity  is  set  accorfling  to 

|u(k)l2  ~  k^expi-2k^/kl)  (17) 

where  k^,  is  the  wave  number  at  which  the  spectrum  peaks;  the  phase  of  the 
u(k)  are  initially  random.  In  setting  the  initial  conditions  on  the  velocity  it  is 
useful  to  decompo.se  it  into  ‘incompressible’  and  ‘compressible’  parts;  in  terms 
of  Fourier  coefficients,  the  decomposition  is  easily  effected: 

u(k)  =  u’(k)-|-u^(k)  (18) 

where  the  solenoidal  {i.e,  incompressible)  part  u’  and  compressible  part  u''  are, 

u'’(k)  =  (I-kk)-u(k) 

u''(k)  =  kkn(k)  (19) 

and  where  k  is  the  unit  vector  in  the  direction  of  k.  The  k  =  0  component  of 
u’  corresponds  to  the  mean  How  velocity. 

For  the  moment  assume  that  we  have  an  initially  incompressible  flow.  Then 
ii'"  =  0  and  ft  —  1  at  I  =  0;  upon  taking  the  divergence  of  (10)  and  defining 
T  =  1  +T',  we  find  the  fluctuation  T'  initially  satisfies; 

-v^r  =  lOiOjti’u’ 

=  .S’(x)  (20) 

In  terms  of  Fourier  coefficients,  the  temperature  fluctuations  which  are  consis¬ 
tent  wit  h  the  assumption  of  an  initially  incompressible  flow  are  given  by 

T'(k)  =  ^•““.S’(k)  where  k  >  Q  (21) 

Since  7’'(x)  >  —  1  in  order  that  T  >  0,  and  since  T'  ~  then  u''  must 

be  scaled  accordingly.  The  mean  value  of  lii’p  (denoted  by  (|u’p))  is  the  ini¬ 
tial  turbident  Mach  number  squared  M^.  The  maximum  fluctuation  of  |u''p 
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is  larger  than  Mf ,  and  we  expect  that  maxA/j  ~  7“'  in  order  that  the  con¬ 
straint  7’'  >  —  1  is  obeyed.  Here,  when  initial  flow  conditions  with  close  to 
maximal  incompressible  velocities  are  desired,  the  initial  velocity  is  scaled  so 
that  minT'(x)  =  —0.99. 

In  truly  incompressible  initial  conditions,  the  flow  is  subsonic.  However, 
we  are  at  liberty  to  choose  A,  u,  and  cr  arbitrarily  as  initial  conditions  for  the 
Navier-Stokes  ecpiations.  Thus,  supersonic  (max  |u(x)|  >  1)  initial  turbulent 
flow  conditions  can  be  specified.  For  example,  A(k),  u''(k),  »i’(k),  and  o-(k)  can 
be  specified  independently  of  one  another.  As  another  example,  an  initially  in¬ 
compressible  flow  field  can  be  specified  (to  represent  a  local  region  of  turbulence) 
and  a  compressible  flow  field  (corresponding  to  an  approaching  shock  front)  can 
be  added  to  it.  These  alternatives  will,  in  fact,  enable  us  to  investigate  highly 
compressible  homogeneous  turbulence  and  shock-turbulence  interactions. 


7  Numerical  Results 

The  pseudospectral  logarithmic  variable  method  was  implemented  in  both  a  3D 
code  and  a  2D  code.  The  3D  code  was  used  to  simulate  supersonic  isotropic 
turbulence  on  a  64^  grid,  while  the  2D  code  was  used  to  examine  the  passage  of  a 
region  of  turbulent  flow  through  a  shock  on  a  512^  grid.  These  numerical  studies 
will  be  discussed  following  the  definition  of  some  quantitative  flow  measures. 

It  will  be  useful  at  this  point  to  define  a  number  of  quantities  which  measure 
certain  characteristics  of  a  turbulent  flow  and  its  simulation.  First,  the  mean  of 
a  quantity  Q  averaged  over  its  values  at  all  grid  points  Xi  will  be  denoted  by 

w)  =  (“) 

t 


where  u  =  2  for  2D  and  ri  =  3  for  3D.  The  turbulent  Mach  number  Mt,  average 
wavenumber  dissipation  wavenumber  kp,  microscale  Reynolds  number  /?a 
and  compressibility  index  x  are  then 


Mt 

= 

i. 

^ave 

= 

kp 

= 

((e) 

fix 

= 

X 

= 

(23) 

In  the  dimensionle.ss  system 

adopted  here,  1/  =  /i.  Also,  Xave 

=  'Zir/kave  and  the 

local  dissipation  rate  e  is 

£  =  H 

[(V 

xu)2-l-(/?-b4/3)(Vu)-] 

(24) 

Dissipation  clearly  has  a  part  related  to  vorticity  V  x  11  and 

to  dilatation  V  ■  n. 
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7.1  3D  Homogeneous  Supersonic  Turbulence 

In  the  previously  cited  work  [2-7],  initial  RMS  Mach  numbers  were  never  above 
Mt  —  0.8,  and  usually  were  much  less.  One  possible  reason  for  this  limitation 
is  the  occurrence  of  negative  values  of  density  or  temperature  on  the  numerical 
grid,  which  can  lead  to  such  local  instabilities  as  effectively  negative  diffusion. 
Here,  these  possibilities  are  explicitly  avoided  through  the  use  of  logarithmic 
variables.  The  initial  RMS  Mach  number  for  the  test  cases  considered  ranged 
from  Mt  =  1  to  Mt  =  2  and  there  were  no  problems  experienced  in  these  64^ 
simulations  (such  as  conservation  of  energy)  as  long  as  ihe  dissipation  wave 
number  began  and  remained  less  than  ^mar  =32. 

Clonsider  two  cases  for  which  M,  =  1  initially;  in  these  ca.ses,  half  of  the 
physical  volume  contains  flow  with  supersonic  velocity  at  <  =  0.  In  both  cases, 
7  =  1 .4,  /i  =  0.01 , 0  =  0,  and  k  =  0.7  (the  Prandtl  number  is  thus  approximately 
unity);  these  values  correspond  roughly  to  those  of  air.  Also  at  f  =  0,  A  =  <t  =  0 
everywhere  and  the  velocity  satisfied  (17),  with  k^,  =  6. 

The  difference  between  the  two  runs  lay  in  the  value  of  x-  In  mn  3D64A, 
X  =  0.25  while  in  run  3D64B,  x  =  0.75  at  <  =  0.  Both  simulations  ran  until 
almost  f  =  1  at  about  14  cpu-sec/A<;  the  fluctuation  in  total  energy  was  less 
than  0.2%  for  3D64A  and  less  than  0.5%  for  3D64B.  In  Figures  1  and  2,  the  time 
variation  of  the  quantities  Mt  and  x,  and  in  Figures  3  and  4,  the  time  variations 
of  Rx,  have,  and  kp  are  shown  for  runs  3D64A  and  3D64B,  respectively.  For 
these  runs,  kp  <  kmar  =  32,  so  that  both  are  numerically  well  resovled. 

One  difference  between  the  two  runs  manifests  itself  in  the  evolution  of 
enstrophy  ((V  x  u)^)  and  mean  square  divergence  ((V  •  u)^).  These  quantities 
are  are  presented  in  Figures  5  and  6  for  runs  3D64A  and  3D64B,  respectively. 
Since  /^  =  0  for  these  runs,  a  consideration  of  (24),  along  with  Figures  5  and 
6,  shows  that  dissipation  occurs  primarly  due  to  vortical  motion,  rather  than 
to  dilatational  motion  (although  for  a  short  time  after  the  start  of  run  3D64B 
dilatational  motion  is  actually  more  important). 

Another  difference  is  in  the  relation  between  the  solenoidal  and  compressible 
velocity  spectra,  as  shown  in  Figures  7  and  8,  for  runs  3D64A  and  3D64B, 
respectively.  Figures  7  and  8,  which  correspond  to  t  =  0.57  and  I  =  0.56, 
respectively,  show  that  the  dominant  part  of  the  velocity  spectra  at  the  highest 
values  is  the  compre-ssible  one.  However,  Figure  8,  corresponding  to  run  3D64B 
which  had  an  initial  value  of  0  =  0.75,  indicates  that  the  compressible  part  of 
the  velocity  spectra  is  al.so  dominant  over  the  medium  as  well  as  high  /t-values. 


7.2  2D  Shock-Turbulence  Interaction 

In  addition  to  the  homogeneous  .3D  runs  just  described,  a  set  of  2D  runs  on  a 
512^  grid  were  completed.  These  runs,  2D512A,  2D512B,  and  2D512fl,  differed 
initially  only  in  the  value  of  0  assigned  to  each:  10,  .30,  and  100,  respectively; 
for  all  three  runs,  /<  =  0.001,  k  =  0.7,  and  7  =  1.4.  All  the  runs  began 
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with  identical  initial  conditions,  which  consisted  of  a  region  of  incompressible 
turbulence  satisfying  (17)  with  kg  =  9  filling  the  left  half  of  the  2D  grid  and 
a  shock  wave  in  the  right  half  of  the  grid.  Using  the  shock  jump  conditions 
[10,  p.  335],  density,  temperature,  and  velocity  corresponding  to  a  Mach  2 
shock  were  set  inside  the  right  half  of  the  grid,  while  in  the  left  half  of  the 
grid  the  mean  density  was  <  p  >=  1  and  the  mean  velocity  Wcis  <  >=  2 

(the  frame  of  reference  was  such  that  the  shock  front  was  initially  stationary 
and  the  turbulent  region  weis  moving  from  left  to  right  into  it).  At  the  edges 
of  the  turbulent  region,  there  were  transition  regions  of  twenty  grid  points  in 
the  r-direction  (*.e.,  the  streamwise  direction)  in  which  the  interior  turbulent 
velocity  field  went  smoothly  to  zero  (by  a  squared-cosine  taper).  The  turbulent 
temperature  fluctuations  were  determined  by  (21). 

Similarly,  there  was  a  squared-cosine  taper  from  the  free  stream  values  of 
density,  temperature,  and  velocity  to  their  jump  values,  and  back  again,  over 
the  twenty  grid  points  on  the  outside  edges  of  the  shocked  region.  Thus,  we 
begin  on  the  leftmost  edge  of  the  grid  with  freestream  values,  make  a  transition 
into  a  turbulent  region,  make  another  transition  out  of  the  turbulent  region 
back  into  the  freestream  values  and  then  a  transition  into  the  shocked  region, 
followed  by  a  transition  back  to  freestream  values  on  the  rightmost  edge  of  the 
grid.  This  resulted  in  a  periodic  set  of  initial  conditions  which  could  be  treated 
by  a  Fourier  method.  These  initial  conditions  can  be  thought  of  as  1)  a  one¬ 
dimensional  Mach  2  shock  wave,  and  2)  a  localized  region  of  eddy  turbulence 
placed  into  the  freestream  flow  (and  moving  with  it)  just  ahead  of  the  shock 
front. 

Once  evolution  began,  the  front  of  the  shock  steepened  slightly,  and  the 
density  and  temperature  increased  slightly,  leading  to  a  pressure  jump  ratio  of 
5  across  the  shock  front  (it  was  initially  4.5).  Based  on  this  pressure  ratio  the 
Mach  number  of  the  shock  should  have  risen  to  2.1;  the  shock  front,  which  would 
have  remained  stationary  in  the  reference  frame  chosen,  was  observed  to  move 
forward  at  a  relative  velocity  of  0.1,  in  accordance  with  expectations.  During 
the  same  time,  the  back  of  the  initial  shock  pulse  spread  into  a  rarefaction 
wave.  As  an  example  of  the  time  variation  of  R\,  kp,  and  kave,  consider  Figure 
9,  which  pertains  to  run  2D512B,  and  is  similar  for  the  other  2D  runs.  Since 
kp  <  kmar  —  256,  it.  would  appear  that  the  turbulence  is  sufficiently  well 
resolved. 

The  resolution  of  the  shock  front,  in  turn,  gets  better  with  increasing  values 
of  (i.  Consider  Figures  10,  11,  and  12,  where  spanwise  averages  for  vorticity 
and  pressure  are  shown  for  the  three  runs  for  the  same  time  (t  —  0.567).  In 
Figure  10,  which  corresponds  to  =  10,  there  are  some  oscillations  in  the  shock 
front;  these  oscilations  disappear,  however,  in  the  Figures  11  and  12,  which 
correspond  to  /3  =  30  and  fi  =  100,  respectively.  To  visualize  what  happens  to 
the  turbulence  as  it  cro.sses  the  shock  front,  consider  Fiqures  13  and  14,  which 
show  a  small  section  of  the  grid  for  runs  2D512A  and  2D512B  (/i  =  10  and 
p  =  30,  respectively),  both  at  t  =  0.567. 
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In  Fitiures  II  and  12,  the  vorticity  jumps  by  a  factor  of  about  2.!),  which  is 
also  the  case  for  Figure  10,  if  the  initial  overslioot  is  ignored.  'I'his  is  consistent 
with  linear  theory  [16,  17],  The  shape  of  the  shocks  is  also  consistent  with  the 
results  from  previous  simulations  using  ENO  methods  [18], 

8  Conclusion 

In  this  paper,  two  primary  extensions  of  pseudospectral  F'oiirier  niethods  for 
solving  the  compressilile  Navier- .Stokes  equations  have  been  described.  First,  by 
using  a  ‘logarithmic  variable"  formulation  of  the  basic  equations,  it  was  shown 
that  supersonic  homogeneous  turbulence  could  easily  be  simulated.  Second, 
by  using  realistic  values  of  the  ratio  of  bulk  viscosity  to  shear  viscosity,  it  was 
demonstrated  that  shock  structure  could  be  re.solved  and  shock-turbulence  in¬ 
teractions  could  be  examined  by  direct  numerical  simulation. 

The  logarithmic  variable  method  is  not  an  explicitly  con.servative  formula¬ 
tion,  though  it  was  seen  that  the  global  energy  was  very  well  conserved.  In 
addition,  the  Rankine  -Hugoniot  relation  between  pressure  jump  across  a  shock 
and  upstream  Mach  number  was  also  seeti  to  hold.  Although  these  measures 
indicate  that  the  numerical  method  was  accurately  solving  the  equation.;  of  mo¬ 
tion,  an  exhaustive  study  was  not  performed.  Such  a  study  is  more  appropriately 
dotie  in  a  one-dimension  simulation,  rather  than  the  two-  and  three-dimensional 
simulations  presented  liere,  and  will  be  deferred. 

Note  that  shock  cajiture  and  resolution  is  facilitated  by  the  naturally  occur¬ 
ring  bulk  viscosity  tertn.  It  had  been  found  in  pioneering  work  in  the  numerical 
solution  of  How  i)roblems  involving  shocks  [9],  that  the  introduction  of  an  ‘arti¬ 
ficial  viscosity’  was  necessary  to  ensure  numerical  stability  in  solving  the  Euler 
equations.  Other  viable  stabilization  technicpies  for  the  Euler  ecpiations,  includ¬ 
ing  spectral  filtering  and  smoothing,  have  also  been  developed  [19].  In  addition 
to  these  techniques  for  Euler  Hows,  ‘hyperviscosity’  methods  have  been  ajjplied 
in  solving  compressible  Navier-Stokes  Hows,  for  the  purpose  of  increasing  the 
effective  Reynolds  number  of  a  simulation,  as  well  as  ensuring  stability  [20].  In 
comparison  to  these  a/yerdA/ztir  approaches,  the  novelty  here  is,  again,  that  we 
use  a  physically  motivated  approach  involving  the  bulk  viscosity.  It  may  prove 
useful  to  compare  these  different  methods  in  greater  detail,  but  this  is  beyond 
(he  sco|)e  of  the  (iresent  work. 

Although  the  distribution  of  grid  points  is  not  concentrated  at  the  shock 
front,  this  is  not  a  disadvantage.  In  particular,  the  presence  of  turbiihuice  also 
rerpiires  a  sufficient  number  of  grid  points  for  its  resolution;  thus  it  is  advan¬ 
tageous  to  have  an  even  distribution  of  grid  jioints,  since  small  scale  dynamic 
activity  is  occuring  at  essentially  all  parts  of  the  grid.  In  this  w.ay,  the  small 
scale  structure  of  turbulence  and  the  small  scale  structure  of  shock  fronts  are 
treated  eipially. 

In  the  future,  we  ho|)e  to  incorporate  reacting  Hows  within  the  kind  of  sim- 
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Illations  described  here.  This  is  an  important  aiiplication,  since  a  turbulent, 
reacting  flow,  in  the  jiresence  of  shocks,  occurs  in  many  situations  of  interest 
to  the  aerospace  and  astrophysical  communities.  One  specific  area  of  interest 
is  in  the  supersonic,  turbulent  combustion  which  occurs  in  a  scramjet  engine. 
Another  is  the  inflow  through  the  accretion  shock  around  a  black  hole. 
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Figure  1.  RMS  Mach  number  Mt  and  compressibility  index  x  for  run  3D64A. 


Figure  2.  RMS  Mach  number  Mt  and  compressibility  index  x  for  run  3D64B. 
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Figure  3.  Microscale  Reynolds  number  R\,  average  wave  number  kave, 
and  dissipation  wave  number  kp  for  run  3D64A. 


Figure  4.  Microscale  Reynolds  number  Rx,  average  wave  number  kavf, 
and  dissipation  wave  number  kp  for  run  3D64B. 
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Run  3D64A 


Figure  5.  Enstrophy  and  mean  square  divergence  for  run  3D64A. 


Figure  6.  Enstrophy  and  mean  square  divergence  for  run  3D64B. 
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Figure  11.  Spanwise  averages  for  pressure  and  RMS  vorticity  for  run  20^126. 


f’igure  12.  Spanwise  averages  for  pressure  and  RMS  vorticity  for  run  2D512(  . 
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Spanwise  Grid  Points 


Vorticity  Contours  for  Run  2D512A  at  t=0.749 
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Figure  13.  Vorticity  contours  for  a  small  section  of  run  2D5r2A. 


Vorticity  Contours  for  Run  2D512B  at  t=0.567 
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Figure  M.  Vorticity  contours  for  a  small  section  of  run  2Dr)12H. 
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